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Abstract 

An algorithm is presented that incorporates the tensor form of Maxwell's equa- 
tions in a general relativistic electromagnetic particle-in-cell code. The code 
simplifies to Schwartzschild space-time for a non-spinning central mass. The 
particle advance routine uses a fourth-order Runge-Kutta algorithm to inte- 
grate the four- velocity form of Lorentz force. The current density is calculated 
using the curved space-time of the metric. 
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1. Introduction 

The algorithm described here was developed for plasma simulation in an 
environment around a spinning central mass. The versatility of the algorithm 
allows for calculations without spin. Because the algorithm uses a general metric 
explicitly for the description of the space-time, this algorithm can be used as a 
general relativistic particle-in-cell (GRPIC) code. 

The basic equations used for this new code are described in Section 2. In Sec- 
tion 3, the numerical scheme including discretization, field updates, and particle 
mover. The current deposition scheme is described in Section 4. Initialization 
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and stability criteria are discussed in Section 5. Jet formation is described as 
an example application of this new code in Section 6. The concluding remarks 
are discussed in Section 7. 



2. General relativistic particle-in-cell numerical simulation 

Our algorithm is based on a general relativistic formulation of the EMPIC 
algorithm [1]. We incorporated the physical four- vectors (i.e., velocity, current 
and position) along with the electromagnetic field tensor to simulate the plasma 
particle and field dynamics. 

2.1. Formalism for GRPIC 

The equations which control the development of the particles and fields are 
given by the tensor form of the Maxwell and Newton-Lorentz equations [2^ and 
the Kerr metric [5]. 



F^f = ^J'' (1) 

Pal3;'i + Fp-^-a + F^a-fi — (2) 

where J" is the four-current, is the contravariant electromagnetic field 

tensor 

, r"j, is the Christofel symbol of the second kind for the metric, and is 
the four-velocity of the particle. The Latin and Greek indices have the range 
(1,2,3,4). These equations are derived to calculate the electrodynamics of a 
charged particle moving in a curved space-time. The Kerr space-time in Carte- 
sian coordinates is defined using 

2?77 T''^ 

ds^ = -c^dT^ = g^.^dx^'dx" = dx^ + dy^ + dz^ - [cAtf + —k^ (4) 

+ a^z^ 

^ r^{xdx + ydy) + ar{xdy — ydx) + (r^ + a^){zdz + rdt) 
~ (r2 + a2)r ^ ' 

r^ - (i?2 - - a^z' = (6) 
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where = + + , m is the Schwarzschild mass, ma the angular 
momentum about the z axis, and t;* = dx^/dt. It should be noted that in GRPIC 
simulations each particle represents an ensemble of charged particles with the 
finite shape over the grids [H 13 [6]. The underlying physics of the particle 
motion is governed by the tensor form of the Newton-Lorentz equation. This 
form provides the equation for the acceleration of the particle. The acceleration 
is a function of the space-time curvature defined by the metric and the Lorentz 
force due to the electromagnetic field. The local field is described by the Maxwell 
field tensor. The components of the tensor are calculated using the general 
relativistic form of Maxwell's equations. Using these equations the particles are 
moved and the fields and currents are calculated self-consistently. 

2.2. The computational cycle 

At each timestep, the algorithm solves for the fields and the particle motion. 
This cycle is shown in Fig[T] The cycle initially starts with initial conditions on 
the particle positions and velocities along with the field tensor components at the 
grids. The particle parameters ((7,to,x,u) are known at the particle location. 
The field tensor values are known only at discrete points on the computational 
cells which comprise the simulation space. The interaction between the particle 
parameters and the field tensor is implemented by calculating the current density 
on the computational grid . With the knowledge of the currents the components 
of the field tensor are calculated. The particles are then "pushed" using the 
tensor components interpolated from the grid to the particle. 

The algorithm uses the Newton-Lorentz and Maxwell's equations in tensor 
form. 

3. Field Advance 

Before and after the particles are moved, the B components of the field 
tensor are updated in two half-steps. The E components are advanced a full 
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Figure 1: Schematic flowchart for GRPIC simulation 
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timestep after the particle mover and the current calculation as shown in Figure 
[ij The field tensor [2 is given by P-^^ = -F^* = , F^^ = e^''^B\ where time 
component is given by ct — x^. Here, e-''^' is Levi-Civita tensor [5]. 

The electric and magnetic fields are components of the Maxwell field tensor. 
We use the form of the contravariant tensor [2]. The tensor is given in the 
Cartesian coordinates by 
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Using this form allows the quick transformation of the tensor to a covariant or 
mixed form. Each of the components of the tensor are offset in space using the 
Yee lattice [7' configuration. 
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Figure 2: The Maxwell tensor field components on the lattice. The components of the Maxwell 
field tensor are defined on the face and edge of the computational cube. 



An example of the covariant form of the staggered Yee lattice shown in Figure 
[2] is used for the tensor update portion of the code. The same construction is 
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used for the contravariant and mixed field tensors. 

The grid must be modified to represent the fields in the Maxwell field tensor 
form. In the lattice the components are staggered spatially and temporally 
shifted by ^. 

3.1. Magnetic Field Update 

The magnetic shift Fpj{:x.i) Fi3y{iii + ^ + %-), where /3, 7 = 1, 2, 3. The 
magnetic field is updated by a half-timestep. We use the Yee lattice configu- 
ration and Eq. [2] to obtain a general difference equation for the magnetic field 
components. 

Tjnew Tjold / tp i p ^u^^t, , , 

^al3 — ^afi - (-t^Plia + -fja:l3)^- {>)) 

where a, f3 — 1,2,3 and 7 = 4. Einstein notation does not apply. The discretized 
update equation becomes 

'FJ^^i^ + Sic, J + 52c., k + Ssc) ~ F^^{i,j, k) 



{',j,k) = F^p{i,j,k) 



F" + 5ip,3 + 52p,k + 5^p) - F" (z, j, k) 



^ (10) 



AxP 

3.2. Electric Field Update 

The Yee lattice is also used for the electric field update. The electric shift 
is FaiiTii) — > Faii^ii + ^) , where a = 1,2,3,4. The electric field update is 
governed by the Eq. [l] We write the equation using Einstein notation. All 
Greek indices are in the range (1,2,3,4). Using Eq. [ijthe electric field update 
becomes 

f:L = Ff, + (^^ J" - r:;^F^^ - r^„F"^ - Fj) dx^ (11) 

The discretized update equation becomes 

A-jr 

F^Ui,],k) = F^\i,j,k)~ —J^{i,j,k) 

- r^^(*, J, k)Ff{i, J, k) - jz, j, k)F:'^{i, j, k) 

^ FrihJ, k)) - Fr(» - ^1., J - 52., k - 5s.) ] 4 
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3. 3. Integration of the equations of motion 

Because of the complexities of particle motion in curved space-time we use 
the fourth-order Runge-Kutta (RK4) method to integrate the four-velocity and 
position of each particle. This method is fourth-order accurate in time which is 
crucial for maintaining numerical stability and accuracy as the particles enter 
the region of increased curvature near the central mass. As the values of the 
central mass and spin increase the space-time curvature is greater. 

EMPIC codes generally use the leap-frog method of Boris [8] which is second- 
order. The acceleration of the particle is given by 



Tl.u^u^ = ^F^^uf" (13) 

(14) 
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^ = a«(x,u) = ^F^u^ - Yl^u^u^ (15) 
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7* = /("^ , a J = (16) 



the steps for RK4 at timestep n are: 
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274 

<+i = ^KAti+2(m^At2 + u?At3) + u?At4) 
<+i = ^KATi + 2(&^AT2 + a^AT3)+ajAr4) 

The general relativistic modification of the RK4 algorithm allows for the 
time dilation felt by each individual particle. The physics of the particle motion 
is locally determined by the particle position and velocity at a location in the 
four-dimensional space-time manifold. This insures that the total effect on 
the current and field tensor is due to the cumulative effect of the local space- 
time curvature around the particle. During the particle "push" the Christoffel 
symbols {T'^^) are calculated at the current particle position. The field tensor 
components at the cell nodes are interpolated to the particle using a "cloud- 
in-cell" method The method is modified for curved generalized space-time 
metrics. 
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Figure 3: Schematic of weighting for cell(i,j,k). The shaded portion corresponds to Wij^t 
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The volume of each sub-cell is given by the general equations 

3 2 1 

r^f r^f r^f 

VOLUME{xl,x),xl,x),xl,x))^ j / / V-det(5) dx^dx^dx^ (17) 

^a;^ "^^s "^^s 

where g is the Kerr metric. The weights for each node of the cell is given by 

CELL_VOLUME(z,j,fc) = VOLUME{xlxl+,,xlx^^+^,xlxl+,) 

W,j,fe = VOLUME{xl, x^+i, 4, a;2^i, 4, a;3^i)/CELL_V0LUME(i, j, k) 

W^,+i J, fe = VOLUME{x],xl, xl, x^+i, 4, x3^i)/CELL_V0LUME(i, j, fc) 

W,,,+i,fe = VOLUME{xl xl^,,xlxl 4, a;|+i)/CELL_VOLUME(i, j, k) 

W,,j.k+i - F0L;7M^(4, x^+i, 4, a;2^i, a;^)/CELL_VOLUME(i, J, fc) 

W^,+i,j+i,fe = VOLUME{x\,xl, x^xl, x^, a;^+i)/CELL_VOLUME(i, j, fc) 

W^,,j+i,fe+i = VOLUME{xl a;i+i, a;2, 4, x3)/CELL_V0LUME(i, j, k) 

J, fc+i = V^OL[/M£;(a;i, 4, x% x]+^,xl, a;3)/CELL_V0LUME(z, j, k) 

- V0LUME{x],xl,xlx%xl,xl)IG¥AA.Y0UJU¥.{i,3,k) 

These weights can now be applied to calculate the values of the field tensor 
{Fp) at particle, p, using 

Fip^p = fc)* W,j-fe + f^(i + l,j,fc)* + 

F|(i, J + 1, fc) * + J, fc + 1) * VK,j, fe+i + 

F1(i + 1, j + 1, fc) * + ^^/3"(* + 1, + 1) * fc+i + 

j + 1, fc + 1) * M^^j+i,fc+i + + 1, J + 1, fc + 1) * 

4. Current Deposition Algorithm 

The charge conservation algorithm is based on Umeda et al [5]. It is faster 
than Buneman et al [1] and uses a straightforward coding algorithm. The cur- 
rent density and flux through the faces of the cell calculations are modified to 
incorporate the space-time curvature. The current is calculated using the four 
velocity, J = pu — p-fv. This requires that 7 is calculated for each particle's 
position and velocity. 
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The area of the cell face given by the a;' — ar-'-plane is 

AREA{xl,x),xl,x>f) = j ^ / ^^gugjj- gf^ dx'dx^ (18) 

where Xg and x j denote the start and finish locations of the nodes. The weights 
are given by normalizing each subarea with the total area of the face. This 
process is similar to the volume weighting for the particle mover. The value for 
the surface weights are given by 

FACE_AREA(i,j) = AREA{xl,xl_^_^,x^ ,x-^+i) 

Wij = AREA{x'p,xl+i,x^,xj+i)/FACE_AREA{i,j) 

Wi+ij = AREA{xi,xl,xl,x-^^-^)/FAGF._AREA{i,j) 

Wij+i = AREA{x^j„xl+i,x:l,x^)/FAGE_AKEA{i,j) 

Wi+ij+i = AREA{xl,xl„x'j,x^)/FACE_AREA{i,j) 

The current deposition in the ii'-direction at each node is then given by 

J,^j = FLUX'^ * Wij 

J^.+i = FLUX^'^Wij+i 

5. Initialization and Stability Criteria 

Basic criteria should be set for any PIC simulation. They are related to the 
plasma frequency and the Courant condition. These criteria are designed to 
prduce the essence of the plasma dynamics. Some of these criteria are modified 
and interpreted for general relativity. 
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Figure 4: Area schematic for generalized — 
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5.1. Plasma Frequency 

The plasma frequency condition parameterizes the plasma oscillations of the 
system on the lowest order. The plasma frquency is given by 



where Wp is the plasma frequency, no is the particle density, q is the charge and m 
is the mass of the particle. At is the simulation timestep. These conditions help 
control the numerical stablity of a PIC simulation. The parameters determined 
by the initial conditions of the simulation are n, q, and m. The plasma frequency 
criterion is modified by the general relativistic effects of curved space-time. The 
partilcle density in curved space-time is transformed by n 7710- The time 
dilation caused by general relativity transforms At Ato /j = At. The general 
relativistic plasma frequency to lowest order becomes 



cOj,At < 2 



(19) 



(20) 



coAt = ^/^imxfjm = ■y~^^'^{(jpAt) < 2 



(21) 
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Near the central mass, 7 will increase, therefore a given uq and Ato which 
satisfy the plasma frequency criterion in Minkowski space will be an approximate 
stability condition near the central mass. The same arguement applies for the 
Courant condition. These conditions illustrate the difficulty of using classical 
plasma parameters as criteria for the simulation near the event horizon. 

5.2. General Relativity Considerations 

The coUisionless dynamics associated with relativistic outflows can be stud- 
ied using GRPIC. Speciflcally, the algorithm is well suited for modeling plasma 
environments in strong gravitational fields. The difficulty of incorporating gen- 
eral relativity in PIC algorithms is that the scalings are typically vastly different. 

Daniel and Tajima [TD] used the PIC algorithm in a 3 -I- 1 Schwartzchild 
metric. The simulation was a l| dimension PIC simulation around the event 
horizon. The PIC algorithm is well suited for studying the acceleration mech- 
anism around the event horizon. Also particle codes lend themselves to the 
incorporation of future electromagnetic radiation studies. 

General relativity uses time and length scales on the order of the Schwartzchild 
radius, Rg = GM/c^, and time, = Rs/c. PIC codes generally use scales in 
terms of the plasma skin depth, Ae = c/up, and the period, r ~ Because 
of these vastly different scales a one-to-one global simulation is not feasible at 
this time. However, important physics can be understood by scaling the physics 
to dimensionless units. This allows the exploration of the dynamics of the sys- 
tem. For the GRPIC algorithm we use Stoney scalings [HJ to normalize the 
physical quantities in our simulation space. We set the standards of the system 
to be Ax^ ~ Ax'^ — Ax^ ^ At = 1, c is set to satisfy the Courant condition, 
and the unit charge, q, is used to normalize the charge dimensions. These values 
give the dimensionless gravitational constant and Schwartzchild radius as 




(22) 




(23) 



where Mc is the mass of the central object. 
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6. Application 

The initial disk geometry of our simulation consists of a free falling corona 
and a Keplerian disk (see Figure [5^), which is similar to Nishikawa et al [TSIIT^ . 
The central mass is located at the center of the computational space and is co- 
rotating with a = .95. The central mass has a value of Mc = 300. The Keplerian 
disk is located at r > = 3rs | cos 6*1 < 5, where (5 = 1/8 is the disk radius. 
In this region the particle number is 100 times that of the corona. There are 
0.8 million disk particles. The initial orbital velocity of the particles in the disk 
is V(j, — V}^ = c-y/rs/(2r), where uk is the Keplerian velocity, where r is defined 
by Eq. [6] There are no disk particles initially at r < r^. The initial magnetic 
field is taken to be uniform in the z direction. The magnitude of the field is 30 
field simulation units. This field component is the contravariant z component 
of the field. The charge-to-mass ratio of the particles is 10^'^. Our simulation 
domain has a mesh size of 64 x 64 x 128. The computational space is scaled using 
the nodes. The grid is uniform and normalized to Aa;* = 1, where i — 1,2,3. 
The time scaling is found by using the Courant condition [5j. The speed of 
light, c, is set to 0.05 which ensures that the particles do not move faster than 
the fields. This condition also controls the growth rate of the fields and helps 
prevent nonphysical growth rates during the simulation. 

The jet has a structure which forms spirals around the z-axis. The faster 
particles are at the center of the jet. The number of highly relativistic particles is 
much lower than the slower moving particles within the jet. The initial particles 
forming the jet are moving at high velocity through the infalling corona. The 
jet forms with fast spiraling particles around the central axis. The jet consists 
of a distribution of electrons and positrons. 

Although the corona is falling into the central mass, slower particles are 
being swept up by the jet. The slower particles continue to move with the jet. 
The faster particles are at the center of the jet. The number of highly relativistic 
particles is much lower than the slower moving particles within the jet. 
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(a) Initial disk and (b) Beginning jet (c) Fast spiraling jet (d) Development of 
corona formation slow moving outside 

particles surround- 
ing faster particles 
at the core 



Figure 5: Jet formations are shown in the four panels at the times t = 0, 100, 300, and 500. 
The jets are bipolar. Particle pairs are moving through the jet at different velocities. The 
jet has a structure which forms spirals around the z (central) axis, (a) The disk and corona 
are initialized for the simulation, (b) The initial particles are forming the jet. (c) The jet is 
forming with fast spiraling particles around the central axis, (d) The slower particles continue 
to move with the jet. The faster particles are at the center of the jet. 

7. Conclusion 

The algorithm presented here solves a GRPIC algorithm. It incorporates a 
generalized metric in the numerical calculation of tensor form of Maxwell's equa- 
tions and the Newton-Lorentz equation. If sufficient care is taken to satisfy the 
numerical stability, then this can be a very useful algorithm for using particle- 
in-cell codes for the simulation of astrophysical regimes where the gravity due 
to a central mass plays a role in the dynamics of the system. 
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